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ABSTRACT 

Context. The observation of various ices in cold molecular clouds, the existence of ubiquitous substellar, cold lU globules in planetary 
nebulae and supernova remnants, or the mere existence of comets suggest that the physics of very cold interstellar gas might be much 
richer than usually envisioned. At the extreme of low temperatures (< 10 K), H 2 itself is subject to a phase transition crossing the 
entire cosmic gas density scale. 

Aims. This well-known, laboratory-based fact motivates us to study the ideal case of a cold neutral gaseous medium in interstellar 
conditions for which the bulk of the mass, instead of trace elements, is subject to a gas-liquid or gas-solid phase transition. 

Methods. On the one hand, the equilibrium of general non-ideal fluids is studied using the virial theorem and linear stability analysis. 
On the other hand, the non-linear dynamics is studied using computer simulations to characterize the expected formation of solid 
bodies analogous to comets. The simulations are run with a state-of-the-art molecular dynamics code (LAMMPS) using the Lennard- 
Jones inter-molecular potential. The long-range gravitational forces can be taken into account together with short-range molecular 
forces with finite limited computational resources, using super-molecules, provided the right scaling is followed. 

Results. The concept of super-molecule, where the phase transition conditions are preserved by the proper choice of the particle 
parameters, is tested with computer simulations, allowing us to correctly satisfy the Jeans instability criterion for one-phase fluids. 
The simulations show that fluids presenting a phase transition are gravitationally unstable as well, independent of the strength of 
the gravitational potential, producing two distinct kinds of substellar bodies, those dominated by gravity (“planetoids”) and those 
dominated by molecular attractive force (“comets”). 

Conclusions. Observations, formal analysis, and computer simulations suggest the possibility of the formation of substellar TF clumps 
in cold molecular clouds due to the combination of phase transition and gravity. Fluids presenting a phase transition are gravitationally 
unstable, independent of the strength of the gravitational potential. Arbitrarily small H 2 clumps may form even at relatively high 
temperatures up to 400 - 600 K, according to virial analysis. The combination of phase transition and gravity may be relevant for a 
wider range of astrophysical situations, such as proto-planetary disks. 

Key words. Instabilities - ISM: clouds - ISM: kinematics and dynamics - ISM: molecules - Methods: analytical - Methods: numer¬ 
ical 


1. Introduction 

Typically around 50% or more of the gas in spiral galaxies con¬ 
sists of H 2 , inferred indirectly by CO or dust emission. Since 
the discovery of dark molecular hydrogen ( Grenier et al.|2 005 1 
Langer et al. 2010| |Planck Collaboration et alTj 2011 1 [Paradis 
et al.||2012| >, the estimated quantity of H 2 in the Milky Way 
has essentially been doubled, effectively revealing the nature of 
some of the dark baryons. It is commonly admitted that the CO- 
related molecular hydrogen is present in relatively dense regions 
of the interstellar medium, molecular clouds with number den¬ 
sity > 10 10 m 5 and temperatures of 7-30 K (Draine|2011 


Even though PE is by far the most abundant molecule (~ 
90%), molecular clouds are mainly detected by CO emissions 
because of all the difficulties in detecting cold Hi (Bolatto et al. 
201 3[ >. For example, H 2 only starts emitting at temperatures 
> 512 K. See |Combes & Pfennigerj(T997| i for a review of several 
possible methods for detecting cold H 2 . Because of these detec¬ 
tion difficulties, the real quantity of H 2 (as well as He, which 
shares similar properties of discreteness with H 2 ) in molecular 
clouds is still rather unknown, especially when the gas tempera¬ 
ture is below < 8 K down to the cosmic background temperature 
of 2.76 K. 


The condensation properties of H 2 relevant for molecular 
cloud conditions are well known from laboratory data (Air 
|Liquide || 1976) >. The phase diagram (Fig. [TJ shows the domain 
of pressure conventionally attributed to molecular clouds. One 
has to keep in mind, however, that since molecular clouds are 
highly structured (Pfenniger & Combes ]1994| , and commonly 
observed in a state of supersonic turbulence (Elmegreen & Scalo 


2004), large fluctuations in density and temperature must occur. 


The presence of ice in the interstellar medium consisting 
of heavier molecules, such as H 2 O, CO, CO 2 and NH 3 cov¬ 


ering dust grains, is nowadays well documented (Allamandola 
|et al.|1999] l. Figure[2]shows the location of the critical and triple 
points of abundant molecules existing in the ISM. H 2 ice has 


been detected in the absorption band at 2.417 pm (Sandford & 
|Allamandola|| 1993) |Buch & Devlin|[l994| |Dissly et al.|| 1994^ 7 
but the interpretation of this detection is that H 2 is mixed within 
H 20 -rich grains in conditions that are too warm to allow the bulk 
of H 2 to condense ( |Kristensen et al.|2011[ ). 

High-resolution pictures of nearby planetary nebulae or rem¬ 
nants of supernova have shown the presence of substellar frag¬ 
ments ( [Walsh & Meaburn||1993| >. These very cold globules, or 
knots, each of the size of a few tens of AU, are at least as cold in 
the inner parts (~ 10 K) as molecular clouds, but much smaller. 
An important feature is that the apparent column density in these 
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knots increases inwards as long as the resolution allows, or until 
the knot is optically thin (Burkert & O’dell 1998} . If these trends 
extend to the centre, one can expect there would be much higher 
density at colder conditions. It would be ideal to eventually reach 
a regime where Hi could condense in liquid or solid form, espe¬ 
cially because at high column density the medium blocks UV 
radiation and cosmic ray heating. 

At the level of molecular clouds, even though their average 
properties are well separated from the H 2 phase transition, these 
are only static properties that ignore the highly dynamical nature 
of supersonic turbulence observed in the interstellar medium, 
where fluctuations must be large. Lower temperatures can be 
reached with fast adiabatic decompression alone. An example 
of fast decompression is displayed in the Boomerang Nebula, 
which reaches a temperature of only 1 K because of a fast ex¬ 
panding wind ( Sahai & Nyman|1997} . Thus, one should expect, 
to be coherent with supersonic turbulence, that regions of expan¬ 
sion and compression must be common in the ISM. 

The fragmentation of self-gravitating gas leading to collapse 
has been studied for over a century. The molecular cloud is nor¬ 
mally represented as a self-gravitating, one-phase fluid (i.e. pure 
gas), which is governed by the balance between gravity and gas 
pressure. The gravity of this fluid is directly proportional to its 
density, and the pressure to the adiabatic density derivative of 
the pressure, ( dP/dp) s . It has been shown by Jeans] (1902] ) that 
if a perturbation with a long enough wavelength is introduced 
into the system, this perturbation will grow exponentially. With 
growing perturbation, the pressure of the fluid will not be able 
to withstand its gravity anymore, which will ultimately lead to 
gravitational collapse. 

The approximation of molecular clouds as a one-phase fluid 
is valid in many cases, but when considering very cold, high- 
density regions, it is an over-simplification. In these condition, 
H 2 will start creating ice, which has to be taken into account 
( Walker}2013) . The dynamics of this kind of fluid presenting a 
phase transition is different from a one-phase fluid, the most im¬ 
portant difference being a very low value of ( dP/dp)^ (Johnston 


2014) . But, as ( dP/dp) s is crucial for the stability of a fluid, a 
cold high-density fluid presenting a phase transition is therefore 
expected to be very unstable. 

Based on the observation of substellar globules in planetary 
nebulas and the dynamics of fluids presenting a phase transi¬ 
tion, we may expect the presence of small, substellar hL ice frag¬ 
ments in molecular clouds due to the fragmentation of cold high- 
density regions. It is of great astronomical interest to study the 
nature of this substellar fragmentation, as the resulting bound 
objects may be too small to start nuclear fusion and, thanks to 
their very low temperature, are very difficult to detect. Some of 
the baryonic dark matter may actually consist of these fragments 
(Pfenniger & Combes| 1994] (Pfenniger et al.|1994) . 

In this article, we study the physics of a self-gravitating van 
der Waals fluid presenting a phase transition analytically and 
with simulations. In the analytic part (Sect. [2}, the physics of 
a van der Waals fluid and the related Lennard-Jones (thereafter 
LJ) potential are recalled. We calculate the virial theorem tak¬ 
ing both the gravitational and the LJ potential into account. The 
virial analysis helps to characterize different types of fluids. The 
stability of a self-gravitating van der Waals fluid presenting a 
phase transition is then analyzed. 

In Sect. [3| the molecular dynamics simulator LAMMPS 
(Plimpton 1 11995} > is introduced. By proper scaling of physical 
constants the Coulomb force solver is used to calculate the 
gravitational force, and the short rang molecular force is cal¬ 
culated with the LJ force. We introduce the concept of super- 



Fig. 1. H 2 phase diagram (bold line) and He (dotted line) in cold 
and low pressure conditions. 



Fig. 2. Critical points (bullet) and triple points (y) of common 
molecules in the ISM. As for H 2 in Fig. [I] the sublimation curves 
cross interstellar pressure conditions very steeply in pressure a 
few K below their triple point. For example, CO should essen¬ 
tially be frozen below ~ 20 K, so unable to emit rotation lines. 


molecules to enable us to perform simulations with a total mass 
high enough for the fluid to be self-gravitating. 

The simulations performed are discussed in Sect. [4] First, the 
correctness of the super-molecule approach is tested. Second, 
one-phase fluids (i.e. pure gas) are used to test the Jeans cri¬ 
terion. Third, simulations close to a phase transitions are per¬ 
formed, studying the properties of non-gravitating fluids and flu- 
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Fig. 3. van der Waals phase diagram for a fluid with T t = 0.9. 
Dotted line: gas phase; dashed line: phase transition; solid line: 
solid phase. 


Fig. 5. Hi and van der Waals phase diagrams. Solid line: Hi lab¬ 
oratory data; dotted line: van der Waals vapour curve derived 
with the Maxwell construct. 



Fig. 4. van der Waals phase diagram for a fluid with a temper¬ 
ature (from bottom to top) of 7j- = 0.2, 7' r = 0.4, T r = 0.6, 
7j = 0.8, T, = 1.0, and 7j = 1.2. Solid line: van der Waals 
EOS with Maxwell construct, dotted line: original van der Waals 
EOS. 


ids with a gravitational potential above and below the Jeans cri¬ 
terion. 


2. Physics of a fluid presenting a phase transition 

In this section, we describe the physics of a self-gravitating fluid 
presenting a phase transition. 

2.1. van der Waals equation 

The van der Waals equation is a classical equation of state (EOS) 
of a pure fluid presenting a first order phase transition. This EOS 
describes the macroscopic behaviour of a fluid in which at mi¬ 


croscopic level the molecules are strongly repulsive at short dis¬ 
tances and beyond that weakly attractive over a limited range. 
The van der Waals equation (van der Waals|19l0; Johnston 


2014 1 is a modification of the ideal-gas law, taking the finite 
size of molecules and intermolecular interactions into account. 
It links pressure, density, and temperature as follows: 


k B T n 9 

P = —- an 2 

1 — bn 


(1) 


with a [Pam 6 ] and /?[m 3 ] being constants characteristic of the 
fluid. It can also be expressed in a reduced, dimensionless form 
as. 


Pr = 


87V 

XTT 


3 ni 


( 2 ) 


where P r = P/P c , n r = n/n c , and T r = T/T c . The parameters 
P c , T c and n c are the values of the thermodynamic critical point 
(Kondepudi & Prigogine f1998] |Carey[ 1999 1 . 

Figure[3]shows the phase diagram for a fluid with T = 0.97V. 
One can distinguish the gaseous phase (dotted line) and the solid 
phase (solid line); in between (dashed line), the fluid presents a 
phase transition. There are states where ( dP/dp) s < 0, which 
is thermodynamically unstable. The parts of the dashed curve 
where (dP/dp) s > 0 are metastable because a lower entropy state 
is reached when the fluid splits into condensed and gaseous com¬ 
ponents. The fraction of condensed over gaseous phase grows 
from 0 to 1 from left to right. 

When the two phases coexist, the van der Waals EOS is re¬ 
placed by a constant pressure marked by a horizontal line. This 
constant pressure level is determined by the Maxwell “equal 
area” construct (Clerk-Maxwell p~875 i, demanding a total zero 
P ■ v work for an adiabatic cycle between the fully gaseous to the 
fully condensed state. 

Figure [4] shows the phase diagrams for fluids with a temper¬ 
ature from T = 0.27V to T = 1.27V- The dotted line shows the 
original van der Waals EOS whereas the solid line displays the 
modified law using the Maxwell construct. There is no Maxwell 
construct for T > T c as (dP/dp ) s is always > 0. Lekner (1982) 
provides a parametric solution to the van der Waals and Maxwell 
construct, using A s as variable. 

The phase equilibrium pressure is almost identical with the 
coexistence curve of laboratory data for H 2 over a wide range of 
pressure (Fig. [5]). 
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2.1.1. Lennard-Jones forces 



The van der Waals model of phase transition fluid is conve¬ 
nient for a continuum fluid description, but fails to represent all 
the phenomena around the phase transition, which is the reason 
for correcting it with the above mentioned Maxwell construct. 
Even with the Maxwell construct, however, a local thermal equi¬ 
librium is still supposed to hold. In astrophysical contexts, often 
even these assumptions cannot be granted. For instance, in super¬ 
sonic turbulent situations, ubiquitous in the interstellar medium, 
thermal equilibrium cannot be satisfied since mechanical energy 
propagates faster than thermal energy through pressure. Thus 
any method using quantities, like temperature or pressure, im¬ 
plicitly assumes that a local thermal equilibrium is established, 
which makes their use in the supersonic turbulent regime uncer¬ 
tain. 

A much less demanding model of fluid is provided by molec¬ 
ular dynamics, where the simplest molecule interactions close to 
the van der Waals model in equilibrium situations is provided by 
the Lennard-Jones (LJ) intermolecular potential (Jones |1924) >. 
No local or global thermal equilibirum is required. The LJ po¬ 
tential consists of an attractive long-range term and a repulsive 
short-range term (Fig. [6]), 


®u( r ) — 4— (lt' 2 — r /f) , 

m \ ' 


(3) 


with ro- = r/cr. Its minimum value, located at r cr = 2'^ 6 , is -e/m. 


2.2. Virial theorem for a Lennard-Jones gravitating fluid 
2.2.1. Lagrange-Jacobi identity 

The virial theorem (Clausius jl870| describes the statistical equi¬ 
librium of a system of interacting particles or fluid systems, re¬ 
lating the kinetic and potential energies. In the case of a self- 
gravitating LJ fluid, the LJ potential and gravity combine as a 
total potential <D = <t> G + d) LJ . 

The virial theorem for this new potential can be derived us¬ 
ing the Lagrange-Jacobi identity path, by taking the second time 
derivative of the moment of the polar inertia I = jT, , we 
find. 



7max [K] 

Fig. 7. Minimum mass of isothermal equilibrium curves for LL 
homogeneous spheres. 


where the first right-hand term is twice the kinetic energy £km, 
and the second one is the virial term. For a system near a sta¬ 
tistical equilibrium, both sides of this equation should oscillate 
around 0, so the respective time averages should vanish. 

The LJ potential is a sum of two homogeneous functions Q 
®lj = fih.j.a + < I ) i.j,r of degree -6 and -12 respectively, while 
gravity is of degree -1. So we can use Euler’s theorem of homo¬ 
geneous functions to express the virial term as a sum of poten¬ 
tial energies multiplied by minus the homogeneous degree. The 
Lagrange-Jacobi identity becomes 

1 d 2 1 

2 ^2 ~ 12.£pot,LJ,r + 6ZSp 0 t,LJ,a ^pot,G • (5) 

>0 >0 <0 <0 


2.2.2. Homogeneous sphere 

For homogeneous, finite mass spheres at a given temperature, 
the individual terms of the Lagrange-Jacobi identity read, 

2 E km = — M, (6) 


12£'pot ) LJ,r 
ftEpotXJ.a 


pot.G 


m 


12 C r 
-6 c a 



ecr 12 p 4 


M , 


ecr 6 p 2 


M , 


MV 


5/3 


125 


(7) 

( 8 ) 

(9) 


where the constants c r and c a in the LJ terms have been calcu¬ 
lated by straight summation over the 1.5 ■ 10 12 nearest molecules, 
as given in TablejTjfor both the simple cubic (SC) and hexagonal 
close packed (HCP) or face-centred cubic (FCC) lattices. 


Table 1. Lattice constants. 


Lattice 

Cr 

Ca 

SC 

24.8085962 

33.6076959 

HCP/FCC 

48.5275208 

57.8156842 


1 d 2 / 

2d 7 - 


Y mr 2 i + Y m ’ r ‘ ' f i > 

i i 


a By definition, a homogeneous function / of degree k satisfies 
f(Ax) = A k f(x). 
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Fig. 8 . Virial equilibrium of gravitating LJ homogeneous isother¬ 
mal spheres with H 2 molecules arranged on a HCP/FCC lattice. 
The sphere mass as a function of density is plotted for tempera¬ 
tures ranging from r = T/T max = 10 ~ 2 to 10 2 . 


The virial equation, obtained by setting d 2 //df 2 = 0 in the 
Lagrange-Jacobi identity, contains terms proportional to M ex¬ 
cept the gravity term, which is proportional to M 5/2 . Dividing 
the equation by M we can separate M 2/: \ yielding. 


M 2 ' 3 = 


375 
Anp) 


.1/3 


1 

Gm 


k B T +4c r e 




, (10) 


The corresponding density pf at which fragmentation can oc¬ 
cur at arbitrarily small mass is 



which is of the order of the individual molecule density ra/cr 3 . 

The Lagrange-Jacobi identity therefore allows us to predict 
a density pf and a maximum temperature 7 max below which a 
homogeneous sphere made of LJ molecules can find a gravi¬ 
tational equilibrium at an arbitrarily small mass: in a way, this 
provides a simple model, an explanation, for the reason of form¬ 
ing substellar ice clumps out of a cold self-gravitating gas. 
For FL molecules arranged as SC or HCP/FCC lattice, we find 
Pf = 73.8 kg m 3 and 69.2 kg m~ 3 , respectively. These densities 
are of the order or slightly below the solid or liquid Hi density 
0 80 kg nr 3 ). 

At temperatures T > r max the isothermal equilibrium curves 
have a minimum mass when dM/dp\j = 0, which is expressed 
as 


r 3 rr 3 r ll/2 
M min - K Ghn 4 c 5/2 



where t - T / 7 max and 




11/2 

* 0.1686476934. 


(14) 


(15) 


Figure [7] shows the minimum mass as a function of the tem¬ 
perature for HCP/FCC H 2 homogeneous spheres. It is rising very 
steeply at T > 627 K, but rising much slower after ~ 1000 K. 
Interestingly the mass range includes all the masses below ter¬ 
restrial planet masses for temperatures below H 2 dissociation. 


which is indeed positive if the sum of the terms inside the brack¬ 
ets are also positive. However the term proportional to c a is neg¬ 
ative, and when large introduces the possibility that no positive 
solution for A7 2 ' 3 exists. Solving the bracket terms equal to zero 

for (<r 3 p/m) , we have then a quadratic equation for this term, 
which gives the solutions for the densities po for which M van¬ 
ishes. 


( 


3 \ 2 

O' Po l 

m ) 


1 ± 


V>- 4 |¥ 



(ID 


The + and - solutions are real non-negative when the term inside 
the square root above is non-negative. The maximum tempera¬ 
ture at which M can vanish is thus, 

faL max _ c a nof 

e ~ 4c r ( ’ 

Since the critical temperature for a phase transition in the ab¬ 
sence of gravity is about e/fa, we see that T max can be sub¬ 
stantially larger than the critical temperature. With the values 
given in Table 1 for the constants c a and c r , the maximum tem¬ 
perature below which gravity combined with molecular forces 
enhances fragmentation are T max = 414.3 K and 626.8, K for 
the SC and HCP/FCC lattices, respectively, which are 11.38 and 
17.22 larger than e/fa. 


2.2.3. Condensed bodies 

Different condensed and uncondensed bodies can be identified 
using the Lagrange-Jacobi identity, as summarized in Table [2] 


Table 2. Condensed and uncondensed bodies in a LJ fluid using 
the Lagrange-Jacobi identity. 


Dominating Terms 

Name 

2£kin + 62Tpot,LJ,a 

“molecular gas” 

12iip 0 t 5 LJ, r + 6£pot,LJ,a 

“comet” 

1 22l«p 0 t 5 Lj <r + ^pot,G 

“rocky planetoid” 

2^ kin "I” ^pot,G 

“gaseous planetoid” 


Figure [ 8 ] shows some M(p, T) equilibrium curves of grav¬ 
itating homogeneous sphere whose molecules are arranged on 
a HCP/FCC lattice, which is able to represent the most compact 
sphere lattice at high density. At low density, on the left of the di¬ 
agram, the equilibrium is principally fixed by the attractive part 
of molecules and their temperature, and the exact lattice struc¬ 
ture is irrelevant. This is the domain of uncondensed “molecular 
gas”. 

On the right, there is an accumulation curve at approximately 
p * 100 kg m 3 representing purely condensed H 2 bodies mainly 
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balancing gravity with the repulsive part of the molecular inter¬ 
action. The area between the dashed line, connecting the min¬ 
imum of the isotherms with T > 7' max , and the accumulation 
curve is the domain of “rocky planetoids”. 

At temperatures below 7 max = 626.8 K, the equilibrium 
curves plunge to arbitrarily small masses along two regimes: the 
left vertical asymptotes represent the limit gaseous bodies, and 
the right asymptotes the condensed (solid or liquid) bodies called 
“comets”. The area of the “comets” lies between the dotted line, 
connecting the elbows of the isotherms with T < 7' max , and the 
dashed line. 

At high temperature and low density, to the left of both the 
dotted and dashed line and above the isotherm lines, lie the 
“gaseous planetoids”. These bodies’ equilibrium is fixed prin¬ 
cipally by gravity and temperature. 

Of course real astrophysical bodies are not homogeneous, in 
practice in the small mass regime one can expect bodies made 
of a mixture of condensed state in the core surrounded by a less 
dense gaseous atmosphere, which could be calculated by solving 
the hydrostatic equilibrium, see Pfenniger (2004 1 , where some 
models of isothermal bodies made of H 2 containing a solid core 
and a gaseous atmosphere have been discussed. 


2.3. Gravitational instability 

The linearized wave equation for the isentropic density pertur¬ 
bation p per of a self-gravitating fluid of density p is the following 
Peansfl 902] Weinberg|1972| i: 


<3 Pper / 8P \ 2 . „ 

pp ' 


( 16 ) 


Its solution is a superposition of modes of the form 
exp (i(k ■ x - cot)), where co is the frequency and k = 2zr /A the 
wavenumber, and where A is the wavelength. This leads to the 
instability condition (Jeans|1902)>, 


to 


2 = \jp) ^ ~ AnGp < ° ' 


(17) 


When fulfilled the modes are real-exponential, therefore unsta¬ 
ble. The classical Jeans criterion reads, 




(18) 


Therefore the effective gravitational instability is directly depen¬ 
dent on the EOS of the medium. 


2.3.1. Ideal gas 

In most textbooks about Jeans instability only the ideal gas is 
discussed. In the case of an ideal monoatomic gas, the pressure 
derivative term is 


(?) -y~ 

\ °P Is m 


(19) 


with the adiabatic index y = 5/3. This value is always positive, 
which leads to the classical, ideal monoatomic gas Jeans crite¬ 
rion for instability. 


A > Aj = 


rryk B T 

Gpm 


( 20 ) 



Fig. 9. van der Waals EOS, including the Maxwell construct. 


2.3.2. Fluid with phase transition 

When considering a fluid presenting a phase transition the term 
( dP/dp) s may differ from that of an ideal gas. As seen before, in 
the case of the van der Waals EOS, the Maxwell construct must 
be used, which may considerably change the combined stability 
of the fluid (Sect. 2.1 1 . In the presence of a phase transition, the 
EOS gradient is strongly modified, in particular ( dP/dp)T = 0. 
Using 


(?)=-(?) 

\opJ s c v \dp) T 


( 21 ) 


with cp and c v both finite values, we find ( dP/dp) s = 0. 
Therefore, a self-gravitating fluid in a phase transition is also 
automatically gravitationally unstable. 

Figure [9] shows the 3D representation of the EOS including 
the Maxwell construct. Clearly the curved triangular region, the 
phase transition region (on the left) has a very different gradient 
than the almost ideal-gas region at low density and high temper¬ 
ature, or the condensed phase region at high density. It is remark¬ 
able that a temperature drop by one order of magnitude from the 
critical temperature leads to a drop in the critical pressure by 
14 orders of magnitude. For example for PL at 3 K the critical 
pressure is 6.4 ■ 1CT 9 Pa. 


3. Molecular dynamics 


For all simulations, the Large-Scale Atomic/Molecular 
Massively Parallel Simulator (LAMMPS) is used (Plimpton 


19951. The LAMMPS software is a state of the art and widely 


used single and multi-processor code in chemistry, material 
sciences, and related fields. Its abilities to quickly compute 
short and long-range forces and the possibility to use a multi¬ 
timescale integrator make it a suitable tool to perform our 
simulations. 
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3.1. Potential Solver 


The LAMMPS code has a wide range of force fields. For the sim¬ 
ulations in this article we use only the LJ and self-gravitational 
potentials. 

Since the straight calculation cost for exact pairwise interac¬ 
tions is 0(N 2 ), approximate but still accurate methods are im¬ 
plemented that make the calculations possible at a much lower 
cost. A cut-off radius r c is used for the short-range forces. As 
the attractive part of the LJ potential drops with r 6 , it is cal¬ 
culated only for neighbour particles within r c . The neighbours 
for each particle are found using a Verlet neighbour list. This 
list is created with a radius of r n = r c + r s , r s being an extra 
“skin” distance to avoid the recalculation of the Verlet list at ev¬ 
ery time-step. This cut-off method is O (N) and therefore linearly 
scales with the number of particles. 

The gravitational potential drops with r 1 so the long 
range interactions cannot be ignored ._Jt is calculated using the 


Particle 3 -Mesh (P 3 M) method (Hockney & Eastwood 


1981). 


The gravitational potential is split in short-range and long-range 
parts in Fourier space (|Ewald| 1921]): 


& = ( Lsr (l - exp(-k 2 r 2 )) + 6 L r exp(-k 2 r 2 ), 


( 22 ) 


where r s is the splitting distance. The potential <J>s R is calcu¬ 
lated at the same time as the LJ potential, using the same Verlet 
list. The potential <1 >l R is calculated in Fourier space using a fast 
Fourier transform (FFT). 

As LAMMPS is designed for the simulation of chemical sub¬ 
stances not far from terrestrial conditions, and self-gravity is 
too weak to be of any influence, no specific self-gravity mod¬ 
ule is provided. For that reason, a tweak is used by using the 
Coulomb’s potential module for an electric held. 


(f> c - 


CWj_ 

er 


(23) 


where C is the interaction constant, e the dielectric constant, and 
qi j the charges of the molecules. As C is fixed in LAMMPS, 
to correctly calculate the gravitational potential we set the con¬ 
stants as follows: 


e = -1 , 


<h = 



(24) 

(25) 


3.2. Time integration 

The time integration is done using the symplectic leapfrog 
scheme. The drift and kick operators, 

D(At) = x(t + At) — x(t) + Atx(t) , (26) 

K(At) = x(t + At) = x(t) + Atx(t ), (27) 

are applied according to the sequence IXAt) K (yj at 

each elementary time-step. For short-range force, the time- 
step is constant throughout the simulation and set as 10 2 - 
10 3 , which is the typical interaction timescale during nearest- 
neighbour molecular interactions. This is reasonable since 
the repulsive interaction and the finite kinetic energy prevent 
strongly varying accelerations between particles. 

Small changes of individual particle positions can signifi¬ 
cantly change the short-range forces, which is why a very small 
time-step is required. But these small changes have almost no 
influence on long-range forces. For this reason, there is no need 
to calculate the long-range forces at every time-step. 


In regard to short-range calculations, long-range calculations 
are an order of magnitude more costly per step but always need 
the same amount of calculation time: creation of the density-map 
with a fifth order interpolation procedure, a FFT solution of the 
potential, and the same interpolation procedure for finding the 
accelerations. 

With the “reversible reference system propagator algo¬ 
rithm” (rRESPA), a multiple timescale integrator is available for 
LAMMPS that enables us to reduce the number of long-range 
force calculation fTuckerman et al.|1992j ). It enables time inte¬ 
gration in up to four hierarchical levels, but only two are needed 
for the simulations we present. The rRESPA time-integration- 
scheme looks as follows: 





^s R 



^L R 



,(28) 


where n$ R is the number of short-range iterations and set as 10 - 
100, /Vl R is the kick operator using the long-range (FFT) accel¬ 
erations, and Ksr the kick operator using the short-range accel¬ 
erations. 


3.3. Super-molecules 


The maximum number of particles that can be simulated on to¬ 
day’s supercomputers is of the order of 10 9 - 10 10 particles. Even 
using such a huge number of molecules would only result in a 
fluid with a total mass of a few femtograms. It is obvious that 
gravity has no effect on this kind of fluid. For that reason, the 
concept of super-molecules is introduced. This concept is well 
established in galactic dynamics simulation, where super-stars 
weigh typically 10 4 - 10 6 M o , and in cosmology where super- 
WIMPS may weigh as much as ~ 10 67 GeVc 2 particles. Two- 
body relaxation or diffusion due to the low number of super¬ 
particles is negligible, provided the simulation time is not too 
large, depending on the specific problem. Practice has shown 
that this kind of an approximation is valid in galactic dynamics 
for instance if the simulation length is of order 100 dynamical 
times (Binney & Tremain e|2008) l. 

The basic principle of the concept is that each super¬ 
molecule consists of i] molecules, and its mass is therefore. 


'«sm = W«m • 


(29) 


To have the same properties of a super-molecule fluid as for a 
normal fluid, every term of the virial Equ. 0 has to be trans¬ 
formed such that the ratios of the terms are invariant. The kinetic 
energy term reads 

2-Ekin = = iVsM^SM(^g]vi) ■ (30) 

Since the total mass is constant, = -Vsm'«sm, the velocity 

dispersion of molecules and super-molecules is also the same. 


( y SM) _ < y M> • 

From the two LJ terms, Equ. W remain the same if 


, ^-12 
e M o- M 


M 

eucr 6 M 


, ^t-12 

<tSM0- SM 


"‘SM 

£SMC sm 


SM 


(31) 


(32) 


(33) 


Solving these two equations and using Equ. (291, we derive 


^SM - t]e M , 
0 "sm = ^/rjo ~m ■ 


(34) 

(35) 
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Fig. 10. Mean kinetic energy as a function of the maximum dis¬ 
placement in an ice clump. 


The gravitational energy Equ. (|9]> does not change since it does 
not depend on super-molecule properties. 

It is important to ensure that the gravitational force between 
two super-molecules remains small compared to the correspond¬ 
ing LJ force within the short-range force cut-off radius. This is 
assured by setting Fc(r c ) <sc Fu(r c ) with r c being the cut-off 
radius. This leads to the following constraint: 


2 

773 


24e M cr M 


Gin 


2 

M 


( r c,SM 



(36) 


where r Cj sM = (a/csm)- Using a cut-off radius of r c — 4cr and 
hydrogen super-molecules, we find a maximum super-molecule 
mass equal to 5.7- 1CU 6 M®, meaning at least 1.7-10 5 particles are 
required for simulating an Earth-mass, virialized, FU-condensed 
body. 


3.4. Ice clump detection 

The LJ potential has its minimum at r m = 2 l/fi <x. If an ice clump 
is at absolute zero temperature, all bound molecules would 
have this distance from at least one other molecule. But, as the 
molecules in a clump are vibrating, the maximum binding dis¬ 
tance between two molecules is generally larger than r m . 

The simplest clump consists of two molecules. Their mean 
kinetic energy as a function of the maximum displacement can 
be calculated numerically (Fig.[l0|>. 

The mean kinetic energy rises up to the maximum value of 
max = 1.3625, after which it drops again. It can be assumed 
that all stable clumps have a binding distance below this value. 
The distance constraint for two molecules to be bound is there¬ 
fore: 

l cr ,max ■ (37) 

Using ra- max as the threshold distance, LAMMPS provides a list 
of clumps at fixed time intervals. 


4. Simulations 

4.1. Units 

A fluid is defined by the number of super-molecules Asm. the ini¬ 
tial velocity distribution (temperature), the number density, and 
the strength of the gravitational potential. All simulation proper¬ 
ties are molecule-independent, the initial temperature and den¬ 
sity are measured as a factor of the critical values T CI and n a , and 


Table 3. Parameters of the super-molecule test simulations. 


Name 

n/n cr 

T/T CI 

Asm 

SM15 

0.006 

1.5 

4 3 - 200 3 

SM30 

0.006 

3.0 

4 3 - 200 3 

SM60 

0.006 

6.0 

4 3 - 200 3 

SM04 

0.1 

0.4 

25 3 - 200 3 

SM06 

0.1 

0.4 

25 3 - 200 3 

SF" 

0.1 

0.1 

© 

SO 

1 

© 

in 


Notes. All simulations are without gravity and use the basic KDK time 
integration scheme. 

(a) External gravity a = 0.1 (L/r 2 ) 

Table 4. Parameters of the one-phase fluid simulations. 


Name 

n/ricr 

T/T a 

Asm 

7j 

OPO 

io- 2 

1.25 

100 3 

0 

OPGw 

io- 2 

1.25 

100 3 

0.5 

OPGs 

io- 2 

1.25 

50 3 - 160 3 

1.25 


Notes. OPO without gravity and using the KDK time integration 
scheme. All gravitational simulations use the P 3 M gravitational solver 
and the rRESPA time scheme. 

Table 5. Parameters of the phase transition simulations. 


Name 

n/n C! 

T/T a 

Asm 

7j 

OPGs 

10" 2 

1.25 

50 3 - 160 3 

0, 0.5, 1.25 

PT1-1 

10" 1 

0.1 

50 3 , 100 3 

0, 0.5, G = Gopgs 

PT2-1 

io- 1 

0.2 

50 3 , 100 3 

0, 0.5, G = Gqpgs 

PT3-1 

10" 1 

0.3 

50 3 , 100 3 

0, 0.5, G = G OPGs 

PT1-2 

10“ 2 

0.1 

50 3 , 100 3 

0, 0.5, G = G OPGs 

PT2-2 

10" 2 

0.2 

50 3 - 100 3 

0, 0.5, G = Gopgs 

PT3-2 

io- 2 

0.3 

50 3 - 100 3 

0, 0.5, G = Gqpgs 

PT5-2 

10“ 2 

0.5 

50 3 , 100 3 

0, 0.5, G = G OPGs 

PT7-2 

IO" 2 

0.7 

50 3 , 100 3 

0, 0.5, G = Gqpgs 

PT9-2 

io- 2 

0.9 

50 3 , 100 3 

0, 0.5, G = Gopgs 

PT1-3 

io- 3 

0.1 

50 3 , 100 3 

0, 0.5, G = Gqpgs 

PT2-3 

10“ 3 

0.2 

50 3 , 100 3 

0, 0.5, G = Gopgs 

PT3-3 

10“ 3 

0.3 

50 3 , 100 3 

0, 0.5, G = Gopgs 


Notes. All simulations in two different runs: without gravity and pertur¬ 
bation, using basic KDK time integration scheme, and with gravity and 
perturbation, using the P 3 M gravitational solver and the rRESPA time 
scheme. 


the gravitational potential is measured as the factor yj = G/Gj 
of the ideal-gas Jeans gravity Gj, defined as 


G _ 2iyk B T {) 
nm 2 L? 


(38) 


with the box side L = (Asm/h) 1 ^ 3 - It is interesting to note that Gj 
is proportional to temperature and inversely proportional to den¬ 
sity. For all simulations, the time unit is defined as the average 
time of a particle crossing the box 



(39) 


with V 2 = 2>* 3 /A. 

The distance at which the gravitational and LJ forces are 
equal is denoted as 


•RjLJ = FgUgu) - Flj(x'GLj) • 


(40) 
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Fig. 11. LJ potential energy as a function of the number of super¬ 
molecules at t = It. Solid line: SM15"; dashed line: SM30 n ; 
dotted line: SM60 fl . 

(a) See Table 0 


Fig. 13. Z-Coordinate of the centre of mass of “comets” as a func¬ 
tion of time of the SF" simulation. Asm = 50 3 and 64 3 (dotted 
line), 80 3 and 100 3 (dashed line), 128 3 and 160 3 (solid line). 

(a) See Tabled 



Fig. 12. Fraction of bound molecules as a function of the number 
of super-molecules at t = 5r. Solid line: SM04", dashed line: 
SM06". 

(a) See Tabled] 

4.2. Super-molecule concept tests 


Figure 12 shows the fraction of bound molecules with re¬ 
spect to the number of super-molecules. All simulations reach 
the same percentage, independent of the number of super¬ 
molecules. 


4.2.3. Precipitation in an external field 

An external acceleration is applied in the —z direction to a fluid 
within a box with reflecting boundary conditions. The fluid 
is very cold and rather dense (see Table [3J and is therefore 
forming ice grains, or “comets”. Because of acceleration and 
Archimedes’s principle, the “comets” fall to the bottom and stay 
there, similar to snowfall on Earth. 

The purpose of this simulation is to test the timescaling of 
the precipitation as function of the number of particles. Figure 
[13] shows the /-coordinate of the centre of mass of “comets” as 
a function of time. The curves are very similar, with a slight 
over-estimation of the collapse time for the smallest simula¬ 
tion. Therefore the precipitation phenomenon timescale is not 
strongly dependent on the mass resolution. 


4.2.1. Potential energy 

A fluid in a cubic box with periodic boundary conditions and 
initial constant density of n — 0.006 n cr is simulated with 
LAMMPS. Three different initial Maxwellian velocity distribu¬ 
tions with temperature T =1.5, 3.0 and 6.0 T cr , and a number 
of particles from Asm = 30 3 - 200 3 are used. Table [3] shows the 
parameters of the different simulations. 

Figure[TT]shows the LJ potential energy of those three fluids 
with respect to the number of super-molecules. One can see fluc¬ 
tuations in the simulations with low Asm- but the result becomes 
reasonable when Asm > 10 4 and a satisfactory convergence is 
obtained if Asm > 10 5 . Therefore all the subsequent simulations 
are performed with Asm > 10 5 . 


4.3. Fluid with perturbation 

To illustrate the behaviour of a homogeneous fluid perturbed by 
a plane sinusoidal wave, a velocity perturbation in the x direc¬ 
tion is introduced with wavelength A equal to the box side L. 
Starting with a Maxwellian distribution v, of velocities for the 
homogeneous unperturbed case, the x-component of the veloc¬ 
ities is perturbed in such a way as to conserve energy. The per¬ 
turbed x-velocity component for each particle i reads, 

<4 = a [u A! - +/3V sin(wx,)] , (41) 

with u) - 2n/L. The correction factor a < 1 is used to keep the 
same kinetic energy in the x direction, with 


4.2.2. Cluster percentage 

High density, low temperature fluids are simulated during 5r 
at various number of super-molecules. These fluids will form 
“comets” and allow us to compare the final percentage of bound 
molecules. Table[3]shows the parameters of the different simula¬ 
tions. 


a 2 - 


£/ <4 


Zi [Vxi +/3V sin(wx;)] 


(42) 


[5 determines the strength of the perturbation, but is supposed to 
be small enough to be in the linear regime of perturbations. In 
the present simulations, fi - 0.01. 
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Fig. 14. Temperature of the one-phase fluid simulations OPO" 
(straight solid line) and OPGs' 1 as a function of time. Asm = 50 3 
and 64 3 (dotted line), 80 3 and 100 3 (dashed line), 128 3 and 160 3 
(solid line). 

(fl) See Table 0 



Fig. 15. Fraction of bound molecules in “comets” of the one- 
phase fluids OPO (straight solid line) and OPGs as a function of 
time. Asm = 50 3 and 64 3 (dotted line), 80 3 and 100 3 (dashed 
line), 128 3 and 160 3 (solid line). 

4.4. One-phase fluid 

To study the reaction of a pure one-phase ideal-gas fluid far 
from the phase transition but with the above plane-wave per¬ 
turbation, the initial temperature is taken above the critical value 
(To = 1.25 T CI ) and the initial density well below the critical 
value (no = 10 2 n cr ). Three different cases are studied: with¬ 
out gravity, weakly self-gravitating below the classical ideal-gas 
Jeans criterion with jj = 0.5, and sufficiently self-gravitating 
above the ideal-gas Jeans criterion with yj = 1.25. The simula¬ 
tions are performed with Asm ranging from 50 3 to 160 3 . Table 0 
shows the parameters of the different simulations. 

4.4.1. Time evolution 

As the simulations conserve total energy, the formation of 
“comets”, implying a decrease of potential energy, can be mea¬ 
sured by the mean temperature (i.e., kinetic energy) of the sys¬ 
tem. This is equivalent to the release of latent heat when a 
gaseous fluid condenses. 

Figures [14] and [15] display the temperature and the “comet” 
percentage evolutions of the one-phase fluid simulations: OPO 


Fig. 16. Temperature of the one-phase fluid simulations OPGs 
as a function of time. Asm = 100 3 with four different random 
seeds. 



Fig. 17. Density at x = (0.5 ± 0.05)L of the one-phase fluids OPO 
(straight solid line) and OPGs as a function of time. Bold line: 
analytic solution for ideal gas. Asm = 50 3 and 64 3 (dotted line), 
80 3 and 100 3 (dashed line), 128 3 and 160 3 (solid line). 


without gravity and Asm = 100 3 super-molecules and the suffi¬ 
ciently self-gravitating fluid OPGs with Asm ranging from 50 3 
to 160 3 . The weakly self-gravitating fluid of simulation OPGw 
is not shown since it is identical to the no gravity case OPO. 

For the fluid without gravity OPO and the weakly self- 
gravitating fluid OPGw, no particular effect is observed. The per¬ 
centage of “comets” rises to < 2%, which can be attributed to 
the initial fluctuations in the distribution. The small perturbation 
introduced in the x-direction does not change the nature of the 
fluid; its temperature and “comet” percentage remains the same. 

The sufficiently self-gravitating OPGs fluids, on the other 
hand, do change. Their temperatures and cluster percentages are 
rising steeply showing a substantial latent heat release. All sim¬ 
ulations reach the same asymptotic temperature ~ 6T cr , but dif¬ 
ferent initial conditions (i.e., random seeds) yield slightly dif¬ 
ferent behaviours during the collapse. This can also be observed 
in Fig. [16] where four identical parameter simulations but differ¬ 
ent random seeds are run with Asm = 100 3 , leading to a time 
difference of ~ 0.5r for reaching the asymptotic upper value. 

Figure[l7]shows the density increase around the centre of the 
perturbation in the range x = (0.5±0.05)L and compares them to 
the analytic solution for an ideal gas. The simulations follow the 
ideal-gas solution with a slight dispersion up to t « 2.5r, where 
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t [r] 


Fig. 18. Evolution of perturbation density. Straight dotted line: 
analytic solution for ideal gas. Asm = 50 3 and 64 3 (dotted line), 
80 3 and 100 3 (dashed line), 128 3 and 160 3 (solid line). 



Fig. 19. 3D-view of simulation OPGs with /Vsm = 100 3 at t = 4r. 


a density rise is not possible anymore because of the repulsive 
LJ force; the density declines for a short while because of the 
collapse rebound down to a minimum at t « 2.7r. This density 
local minimum corresponds with the break in the temperature 
increase visible in Fig. 14 Subsequently, the temperature growth 
is much slower. 


Figure [TS] shows the evolution of density in the simulations, 
and the predicted growth rate line for an ideal gas subject to 
Jeans’ instability (Sect. |2.3[ >. All curve slopes correspond ini¬ 
tially to the predicted Jeans’ growth rate. The initial steep growth 
and the small time delay are caused by the initial noise in the 
particle distribution. Thus, a certain time is needed for the per¬ 
turbation to overcome the noise. 


6 

5 

4 

3 

2 

1 

0 
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Fig. 20. Distribution of y- and z-coordinates of “planetoid” cen¬ 
tre of mass on a 10x10 grid; all OPGs simulations at t — 2.5r. 
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Fig. 21. Temperature distribution of unbound molecules and 
“comets” as a function of “comet” mass of OPGs with Asm = 
100 3 at t = 4t. 




4.4.2. “Comet” mass distribution 


Figure 33 shows snapshots and Fig. 34 the “comet” mass dis¬ 
tribution at different times. During the temperature rise, the 
“comets” are distributed according to a power law. 


— £c for A/comet ^ If u 


(43) 


d log ( jjj M comet) 

dlog (M comet ) 

and no “planetoid” is formed. The power-law index is very steep 
at the beginning with £ c ~ -10, but quickly decreases as it 
reaches a value of * -2.5 at the end. 

After the temperature increase break at t w 2.7r, the power- 
law distribution remains for small “comets”, but one bigger 
“gaseous planetoid” with over 1% of the total mass forms. It 
is shown in Fig. [19] The spherical body seen at 1 < 3r in the 
snapshots is uncondensed gas, which is why it does not figure 
in the “comet” mass distribution (see Sect. |3.4| for a discussion 
how condensed matter is identified). 

A second power law: 


d log ( comet) 

d log (M comet ) 


— (,'p for comet 3 -Tsm 


(44) 


on the right side of the diagram describes the “comets” and 
“planetoids”. Each body occurs typically only once with our lim¬ 
ited Asm. and = 1. Equation (431 describes the mass distribu¬ 


tion of small bodies where gravity is negligible whereas Equ. 


11 























A. Ftiglistaler and D. Pfenniger: Substellar fragmentation in self-gravitating fluids with a major phase transition 



Fig. 22. Mean temperature over four one-phase fluid simulations 
OPGs with different random seeds as a function of time. Asm = 
50 3 and 64 3 (dotted line), 80 3 and 100 3 (dashed line), 128 3 and 
160 3 (solid line). 



n sm 

Fig. 23. Time and size of the largest “comet” at the appearance 
of the turning point. 


( |44| describes large bodies where gravity is weak, but cannot be 
neglected. 

The first high-density plane is created thanks to the plane 
perturbation parallel to the yz -plane at x — 0.5. One can ob¬ 
serve in Fig. [33] at t = 2r that two filaments form, along the y- 
and z- axes, connecting the “planetoid” with itself thanks to the 
periodic boundary conditions. They are aligned with the mesh, 
but not primarily due to grid effects, but because these lines are 
the shortest distance between the “planetoid” and its periodic 
images. As can be seen in Fig. [39] diagonal filaments are also 
possible using the second shortest distance. 

The position of the “planetoid” is, thanks to the plane-wave 
perturbation, situated in the x = 0.5 L plane, the y and z positions 
are random and depend on the initial condition and vary with 


every simulation as shown in Fig. 20 


Figure[2T]shows the temperature distribution in the “comets” 
as function of mass. As the number of large “comets” with the 
same number of super-molecules is small or even one, no error 
bars can be seen for the larger “comets” and the “gaseous plane¬ 
toid”. 


The broad Maxwellian velocity distribution is seen for the 
unbound molecules (M comet = 10 6 ). The small “comets” (2 or 
more particles) clearly have a lower temperature and dispersion 
than the average temperature of the system. This is due to the 
fact that in order for two super-molecules to cluster together, a 
third one is needed, and this takes away some of their kinetic 
energy. The larger a “comet” is, the closer its temperature is to 
the system average temperature. On longer formation time the 
“planetoids” temperature tends towards the system average tem¬ 
perature, as expected. 


termolecular interactions, the simulations are therefore overem¬ 
phasizing the gravitational effect. 

Figure [35] shows the “comet” mass distribution of all OPGs 
simulations. For small “comets”, the percentage for a same num¬ 
ber of super-molecules is the same for all Asm- For example, at 
t — 5 t, all simulations have a fraction of ~ 10 1 for comets con¬ 
sisting of two super-molecules. These small clumps should be 
called multimers as mentioned below. But, with increasing Asm, 
the mass of a comet consisting of two super-molecules dimin¬ 
ishes since AT>sm = 2/Asm- 

On the other hand, for large “comets” and the “planetoid”, 
the mass is invariant to Asm- For example, the “planetoid” at 
t — 5t has the same mass of « 5 ■ 10~ 2 M tot for all Asm- 


The fact that the mass distributions are shifted to the left 

Only the two 


for larger Asm is misleading, as shown in Fig. 36 


largest simulations are shown, the simulation wit 


1 Asm - 128' 


is shown normally, but for the simulation with Asm = 160 3 is 
downscaled to 160 3 /2 ~ 128 3 , always two “comet” sizes are 
added together. As can be seen, when comparing the two simu¬ 
lation using the same scaling, the two mass distributions are in 
fact very similar. 

It is interesting to look at the turning point, defined as the 
point where the mass sum of the biggest “comets” is larger than 
the mass sum of smaller comets. This point is interesting as it 
is the first indicator of the creation of a aggregate of molecules 
that start to be influenced by gravity. Figure [23] shows the time 
and the mass of the turning point. One can see than the time of 
appearance is independent of Asm- The “comet” mass, on the 
other hand, clearly follows a power law with: 


Tfcomet 

M tot 


10 A! 


SM 


(45) 


4.4.3. Scaling 

Figure [22] shows the mean temperature value over four simu¬ 
lations with different random seeds. One can see that the two 
smallest simulations with Asm = 50 3 and 64 3 are starting to col¬ 
lapse earlier than the other simulations, which are very similar. 
This can be explained by the fact that these simulations have 
very low xqlj values (Equ. |40|), 3.65 and 4.03 r, r . The first value 
is in fact lower than the cut-off radius of 4 r, r and fails to satisfy 
Equ. ( [36| , the other one just slightly above it. In these simula¬ 
tions, the gravitational forces are strong even in short-range in- 


with ~ -1- Since Asm = (Wtot/TfsM- at turning point the num¬ 
ber of super-molecules is always of order of 10. In other words, 
in critical conditions where the attractive molecular force adds 
up to gravity, the tendency to make a larger condensed body 
mass fraction already starts at about ten molecules. 

4.4.4. Extrapolation to physical scale 

Having studied the scaling behaviour of the simulations, we 
can attempt to extrapolate some properties for a fluid consist¬ 
ing of real molecules instead of super-molecules. Considering 
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Fig. 24. Extrapolation of the “comet” mass distribution for a real 
Ht molecular fluid at t = 5 t. 


H2 molecules, the total mass of the OPGs simulations is 1.4 M@, 
the box length equal to 2.3 • ICE 3 AU and r = 1.5 • ICE 2 years. 

Thus, in realistic conditions there would be a total amount 
of 2.5 • 10 51 H2 molecules, the number of molecules per super¬ 
molecule rj varying from 2.0 • 10 46 for Asm = 50 3 to 3.2 • 10 44 
for Asm = 160 3 . So, even for the largest simulations, there is a 
difference in number of particles of more than 40 orders of mag¬ 
nitude with real situations. For that reason, the extrapolated data 
has to be treated with caution. As long as no other physics en¬ 
ters in the interval of scales, this is however not extraordinary in 
regards to common practices in simulation works, such as simu¬ 
lating stars with SPH particles. While a single SPH particle can 
at best represent a mass fraction of 10 6 - 10 1,1 of the total, an 
SPH particle is supposed to behave as the smallest mass element 
in local thermal equilibrium, say a few 1000 protons over the 


stellar mass: ~ ICE , which lies therefore well over 10 4U times 
the SPH simulation resolving capacity. 

As can be observed in Fig. 


15 the fraction of bound 


molecules does not change when increasing the number of 
particles and will remain unchanged for a fluid consisting of 
molecules. The same is true for the size of the resulting “plane¬ 
toid”, as can be seen in Fig. [35] Therefore, we can assume that 
after 5r, the system will consist of ~ 80% unbound molecules 
and a “planetoid” with a mass of ~ 0.07 M®. 

At t — 5r, there is a fraction of 10 1 of two molecules 
aggregates i.e. dimers with a mass of 6.6 • 10 27 kg. These 
multi-molecules clumps should be called “multimers” instead of 
“comets”. As mentioned above the turning point occurs at ~ 2.4r 
consistently for approximately 10 molecules, so its mass is about 
20 m H = 3.3 • 10- 26 kg. 

Using the above values and the two power laws of Equ. f43] 
—[44]), Fig.|24]shows a schematic diagram showing the mass dis¬ 
tribution at t — 5r. While the mass fraction of the turning point 
multimers is tiny, the power-law distribution rapidly increases 
the mass in larger grain- and comet-sized bodies until they be¬ 
come a planetoid mass. 


4.5. Phase transition 

We study several physical conditions close to a phase transition. 
Table [5] summarizes their properties. Three different densities, 
ICE 1 , 10 2 and 10 3 n cr , and three different temperatures, 0.1, 
are chosen. As the one-phase fluid studied in 


0.2 and 0.3 T cr , 
Sect. 


4.4 


also has a density of 10 1 n cr , three additional temper¬ 
atures, 0.5, 0.7 and 0.9 T cr , were used for this density to make 



Fig. 26. Temperature of non-gravitating fluids as a function of 
time. Dotted line: n = lCE'n cr ; solid line: n = lCE 2 n cr ; dashed 
line: n = lCE 3 « cr . All simulations with Asm = 50 3 . 


a link between the one-phase fluid and the fluids studied in this 
section. 

Fooking at the phase transition diagram (Fig. [25) , all fluids 
with T < O.377 are on the Maxwell line. Those with p < 10 2 p c 
are on the ( dP/dp) s > 0 part of the van der Waals phase di¬ 
agram, which corresponds to metastable states on the verge of 
a deep phase transition. The fluids with p = 10 1 p c are on the 
(dP/dp ) s < 0 part of the van der Waals phase diagram, which 
corresponds to unstable states in a phase transition. The fluids 
with T > 0.57' c are still in the stable regime. 

Three different cases are studied: without gravity, strongly 
self-gravitating above the Jeans criterion, and weakly self- 
gravitating below the Jeans criterion. 


4.5.1. Fluid without gravity 


To study the evolution of fluids presenting a phase transition 
without gravity, the xqlj value does not have to be considered 
and the number of super-molecules is set to Asm = 50 3 . Figure 


26 shows the temperature evolution of some of these fluids. In 
the beginning, the molecules merge into small “comets”, which 
leads to a decrease of the potential energy and therefore a rise 
of the temperature. Once the temperature is high enough, the 
kinetic and potential energy reach a equilibrium and the fluid 
remains stable. 


The number of “comets” formed is dependent on the den¬ 
sity and inversely dependent on the temperature. The higher the 
number of formed “comets", the longer it takes for the fluid to 
reach a stable regime. For example, the very cold fluids PT1-2 
and PT2-2 only reach an asymptotic value around t ~ 100t. 

The “comets” remain at moderate size as can be seen in 
the snapshots (Fig. [37) and in their mass distribution (Fig. [38): 
The number of super-molecules per “comets” is < 10 2 Asm. 
The mass distribution follows the power law of Equ. ( |43) , but 
at ~ 10 4 , the mass distribution rises again, presenting a big 
number of “comets” between 10 4 to 10 2 super-molecules. 
Fig. [27] shows a close-up 3D view of a medium-sized “comet” 
and smaller aggregates. 

Figure[28]shows the temperature distribution as a function of 
the “comet” mass. As in simulation OPGs, the temperature for 
small bodies lies below average, but reaches average value for 
more massive bodies. 


13 













A. Ftiglistaler and D. Pfenniger: Substellar fragmentation in self-gravitating fluids with a major phase transition 



p[pc\ 

Fig. 25. Position of the LJ simulations in a corresponding van der Waals phase diagram. Solid line: Maxwell construct for the van 
der Waals EOS; dotted line: van der Waals EOS. Star: fluid presenting a phase transition; bullet: one-phase fluid. 


4.5.2. Self-gravitating fluid above Jeans criterion 


All sufficiently self-gravitating simulations use the same gravita¬ 
tional potential, G = 1.25 Gj (T = 1.25 T cr , n = 0.01 n cr ). As the 
temperatures and densities differ, Gj and therefore j\ is different 
for each simulation, but this parameter is always > 1. 


As jj > 1 for all simulated fluids, they are unstable if per¬ 
turbed. The lower the initial temperature of a fluid, the bigger 
is yj, which translates to a faster exponential growth. This can 
be seen in Fig. 29 where the fluids with a density of ICE 2 « cr 
are starting to collapse one after another, from the lowest to the 
highest temperature. 


The evolution of the fluids is identical to the fluid OPGs: an 
exponential rise of the temperature until it reaches a tempera¬ 
ture break, after which the temperature rises only slowly and a 
single “gaseous planetoid” is formed. Figures [39| and |40| show 
snapshots and “comet” mass distribution, again very similar to 
the fluid OPGs. First, there is a formation of small “comets” fol¬ 
lowing the power law of Equ. and then the “comets” are 
massive enough to attract each other and merge into one big 
spherical “planetoid” (t > It), following the power law of Equ. 

The power-law index also evolves in a similar fashion as 
in simulation OPGs, starting out very steep and reaching a final 
value of > -2. 
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Fig. 27. Close-up 3D-view of a medium sized “comet” and 
smaller aggregates. 
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Fig. 28. Temperature distribution of unbound molecules and 
“comets” as a function of “comet” size of PT2-2" without grav¬ 
ity and Asm = 50 3 at t = 100r. 

(a > See Table H] 



4.5.3. Self-gravitating fluid below Jeans criterion 


To compare weakly self-gravitating fluids in a phase transition 
(on the Maxwell line) with out-of-phase-transition fluids (below 
the Maxwell line), yj is set to 0.5 for all fluids. This means that 
the individual G value is different for each simulation. 

As shown in Sect. |4.5.1 the perturbation of the weakly self- 
gravitating one-phase fluid OPGw does not create any effect 
if yj < 1. Figure [30] shows the temperature as a function of 
time of the self-gravitating fluids PT5-2, PT3-2, and PT2-2, and 
compares them to the non-gravitating fluids. The out-of-phase- 
transition fluid PT5-2 remains unchanged and does not differ 
from the non-gravitating fluid, identical to OPGw. 

The temperatures of the fluids in a phase transition PT3-2 
and PT2-2 are exponentially growing and differ clearly from the 



Fig. 29. Temperature of sufficiently self-gravitating fluids as a 
function of time. Dotted line: OPGs fl ; dashed line: PT1-3"; solid 
line (from left to right): PT3-2", PT5-2", PT7-2", PT9-2 a ; dash- 
dotted line: PT1 -1All simulations with Asm = 100 3 . 

(a) See Table E| 



Fig. 30. Temperature of weakly self-gravitating fluids as a func¬ 
tion of time. Dotted line: PT2-2 and PT3-2, with Asm = 50 3 and 
without gravity; dashed line: PT2-2, yj = 0.5 and Asm = 80 3 
and 100 3 ; solid line PT3-2, yj = 0.5 and Asm = 80 3 and 100 3 ; 
dash-dotted line PT5-2, yj = 0.5 and Asm = 100 3 . 


non-gravitating ones. In the beginning, the LJ forces dominate 
the gravitational forces and the fluids behave like non-gravitating 
fluids. PT2-2 reacts much faster than PT3-2, thanks to its low ini¬ 
tial temperature and fluctuations due to the Maxwellian velocity 
distribution,and forms many small- and medium-sized “comets” 
(see Sect. |4.5.T| . Once these “comets” are formed, their mass is 
high enough to attract each other with gravity, forming a single 
“rocky planetoid”. 

This can be seen in the snapshots (Fig.[4l| and “comet” mass 
distribution (Fig. [42| . During the first 15t, the fluid is identical 
to the low-temperature, non-gravitating fluid (compare to Fig. 
[38| . Only then does gravity start to play a role, one can see the 
formation of a “planetoid” and how it swallows the small- and 
medium-sized “comets” during its growth. 

With its higher temperature, PT3-2 forms almost no small- 
or medium-sized “comets”, which can be seen in the snapshots 
(Fig. [43| and “comet” mass distribution (Fig. |44[ > for the first 
~ 50t. Only after 50r does gravity show its effect with the ap¬ 
pearance of a medium-sized “comet”, which is too massive to fit 
into the power law (-) (see Fig. [43] t = 50t). From then on, this 
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Fig. 31. Temperature distribution of unbound molecules and 
“comets” as a function of “comet” size of PT3-2 with jj = 0.5 
and Nsm = 100 3 at t = 200r. 




R [a] 


Fig. 32. Density of “planetoid” as a function of the radius. 
Dashed line: OPGs with Asm = 100 3 at t - 4r, the “gaseous 
planetoid” consists of 41559 super-molecules (0.042 M tot ). 
Dash-dotted line: PT3-2, with yj = 0.5 and (Vsm = 100 3 at t - 
200t; the “rocky planetoid’ consists of 74208 super-molecules 
(0.074 M to t). Solid lines: hydrogen, higher value: solid; lower 
value: liquid. Dotted line: pf of Lagrange-Jacobi identity, higher 
value: SC; lower value: HCP/FCC. 


“comet” attracts super-molecules and thus grows in size until at 
t ~ 125r, a single big “rocky planetoid” is formed. 

Figure [3T] shows the temperature distribution as a function 
of the “comet” mass of PT3-2. The temperature of the “comets” 
containing few super-molecules lies below the average, but the 
“planetoids” temperature is the same as the average. 

PT2-2, PT3-2 and PT5-2 all scale very well. The two fluids 
in a phase transition, PT2-2 and PT3-2, have an almost identical 
timescale for the exponential growth for Asm = 80 3 and 100 3 . 
The slight difference is due to the initial random seed. As xglj 
depends on yj and the temperature, both much lower than for 
OPGs, its value is clearly above the cut-off radius, which is 7.63 
and 8.34 for PT2-2 and 7.03 and 7.69 for PT3-2. 

PT5-2 shows no difference for any number of super¬ 
molecules, with Nsm ranging from 50 3 to 100 3 . 


4.6. “Planetoid” densities 

Figure [32] shows the densities of the “planetoids” as a function 
the radius of the simulations OPGs and PT3-2 with y = 0.5yj, 
and compares them to hydrogen laboratory data and values 
found in the Lagrange-Jacobi identity (Sect. [272]). 

The density of the OPGs “gaseous planetoid” is below that 
of a conventional solid or liquid. As their temperature is above 
the critical value (T > 6 T CI ), super-molecules are not able to 
condense without the aid of gravity, and because of their high 
kinetic energy, are vibrating at much higher amplitudes. This 
leads to more space between the super-molecules and therefore 
a lower density compared to super-molecules below the critical 
temperature. 

The density of the “gaseous planetoid” drops with increas¬ 
ing radius and hence does not qualify as a homogeneous sphere, 
but its average density is below the pf value of the Lagrange- 
Jacobi identity (Equ. [L3j), in accordance with Tab. 0 and Fig. 0 
for “gaseous planetoids”. 

The density of the PT3-2 “rocky planetoid” lies between the 
liquid and solid phase. There is no continuous density drop as 
for the OPGs simulation, the density remains stable up to almost 
the outer radius and the body can be approximated as a homoge¬ 
neous sphere. Its density lies above pf in accordance with Tab.0 
and Fig. [8]for “rocky planetoids”. 


5. Conclusions 

We used analytic methods (Sect. 2) and computer simulations 
(Sects. 3-4) to study substellar fragmentation of fluids present¬ 
ing a phase transition. The motivating astrophysical context are 
molecular clouds where Hi forms the bulk of the gravitating 
mass and is not very far from condensation conditions. 


5. 1. Analytic results 

The study of the virial theorem, using the gravitational and the 
Lennard-Jones potential energies, has shown that there is a max¬ 
imum temperature below which a fluid can fragment at arbitrary 
small masses. This temperature is an order of magnitude above 
the critical temperature. This shows, granted the right circum¬ 
stances, that “comets” can form at a temperatures an order of 
magnitude larger than the critical temperature. In the case of FL, 
this maximum temperature lies in the range of 400 - 600 K, de¬ 
pending on the solid crystalline structure. 

A van der Waals fluid can be in three different states: 
gaseous, solid/liquid, or in a phase transition where the two 
phases coexist. The latter is defined as lying on the line of the 
Maxwell construct, where (dP/dp) s = 0. A fluid is gravita¬ 
tionally unstable if an introduced perturbation has a wavelength 
above a certain value, which depends on (<9 P/dp) s . Since this 
quantity vanishes for a fluid presenting a phase transition, such 
a fluid is also gravitationally unstable. 


5.2. Simulation results 

We performed simulations using the state-of-the-art molecu¬ 
lar dynamics simulator LAMMPS. We used super-molecules to 
simulate gravitational and molecular forces together with a com¬ 
putationally tractable number of particles. 
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5.2.1. Super-molecules 

We tested the super-molecule concept thoroughly. We achieved 
good scaling for non-gravitational effects if the number of super¬ 
molecules is large enough (> 10 5 ). The sticking point is to en¬ 
sure that the molecular LJ forces are dominating the gravitational 
forces in close-range interactions. This is achieved by setting the 
number of super-molecules to be large enough so that the dis¬ 
tance at which the two forces are equal is above the cut-off ra¬ 
dius 4 ctsm- For an H2 fluid this sets a maximum super-molecule 
mass equal to 5.7 • 10 6 M®, limiting the total mass that can be 
simulated by the available computing power. 

In principle, the super-molecule concept should not perfectly 
reproduce the time dependence of diffusive properties, as bigger 
particles introduce faster relaxation. But our experiments using 
different resolutions spanning several orders of magnitude did 
not reveal important modifications in the timings of major col¬ 
lapse and asymptotic evolution state. This enables us to study 
the fragmentation of large bodies, which we call “comets” and 
“planetoids”. 

5.2.2. One-phase fluid 


This study is general enough to be relevant for many astro- 
physical situations. In the case of H2 as main component, it con¬ 
cerns situations where temperature may drop well below 10K, 
such as in cold molecular gas in the outer galactic disks, in ex¬ 
panding winds of planetary nebulae or supernovae, where cold 
substellar mass condensations are known to form, or in proto¬ 
planetary disks. Other molecules besides H2, such as CO, can 
condense even sooner, but as their mass fraction is low they 
should not significantly perturb the gravitational balance, while 
He should not condense at all and should keep a minimal amount 
of gaseous phase |^] The precipitation of formed “comets” and 
“planetoids” in a gravitational field, separating the gaseous and 
condensed phases, is also a fascinating aspect that is as yet im¬ 
possible to capture with traditional hydrodynamical codes, but 
possible with the molecular dynamics approach. We will pursue 
further simulation work, including more specific astrophysical 
applications. 

Acknowledgements. This work is supported by the STARFORM Sinergia 
Project funded by the Swiss National Science Foundation. We thank the 
LAMMPS team for providing a powerful open source tool to the scientific com¬ 
munity. We thank the referee for a thorough reading of the manuscript and con¬ 
structive comments, which substantially improved the paper. 


We applied a plane sinusoidal perturbation into a one-phase fluid 
with a temperature above the critical value. The results repro¬ 
duce the ideal-gas Jeans instability: no collapse is seen for con¬ 
ditions below the Jeans criterion, while an exponential growth of 
the perturbation is observed for conditions above it. 

As a result, the temperature rises, and small- and medium¬ 
sized “comets” form, some of which later merge into one big 
“planetoid”. An interesting observation is the mass distribution 
of these “comets”, which follows a power law for the small- and 
medium-sized ’’comets”, while the “planetoid” and the largest 
“comets” follow a different power law. 


5.2.3. Phase transition fluid 

We simulated fluids with temperatures below the critical value 
and close to the effective phase transition for three different 
cases: without gravity, sufficiently, and weakly self-gravitating 
(above and below the ideal-gas Jeans criterion). 

Because of the Maxwellian velocity distribution fluctua¬ 
tions, the non-gravitating fluids form small- and medium-sized 
“comets” until the potential and kinetic energy reach an equilib¬ 
rium; thereafter they remain statistically stable. In the absence 
of any long-range force, no “planetoid” forms. 

The self-gravitating fluids with a gravitational potential 
above the ideal-gas Jeans criterion do not react differently to a 
plane sinusoidal perturbation from a one-phase fluid: the pertur¬ 
bation grows exponentially and its temperature rises, and small- 
and medium-sized “comets” form, which ultimately merge into 
one “planetoid”. 

The analysis predicts that fluids presenting a phase transition 
are unstable even if the gravitational potential alone is below the 
ideal-gas Jeans threshold. The performed simulations of weakly 
self-gravitating fluids in the phase transition regime did repro¬ 
duce the prediction. Because of the weak nature of the gravita¬ 
tional potential, the timescale to see this reaction is two orders 
of magnitude longer than for sufficiently self-gravitating fluids, 
but, in the end, a big central “planetoid” forms anyway. While 
the out-of-phase transition fluids do not amplify perturbations, 
those in a phase transition indeed amplify them. 
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Fig. 33. Snapshot time sequence of the simulation OPGs with IVsm = 10()\ showing the smooth formation of a “planetoid”. The 
slice size is 0.025 X 1 X 1 L 3 at x = (0.5 ± 0.0125)L. 
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Fig. 37. 0.2 x 1 x 1 L 3 snapshots at x = (0.5 ± 0.1 )L of simulation PT2-2 with N$ M = 50 3 , and without gravity. 
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Fig. 38. “Comet” mass-distribution snapshots of simulation PT2-2 with ;Vsm = 50 3 , and without gravity. 
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Fig. 40. “Comet” mass-distribution snapshots of sufficiently self-gravitating simulation PT1-3 with Nsm — 100 3 . 
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Fig. 42. “Comet” mass distribution of weakly self-gravitating simulation PT2-2 with Nsm = 100 3 . 
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Fig. 43. 0.025 x 1 x 1 L 3 snapshots at x = (0.5 ± 0.0125)L of weakly self-gravitating simulation PT3-2 with Asm - 100 3 . 



-^fcomet/^kftot 

Fig. 44. “Comet” mass distribution of weakly self-gravitating simulation PT3-2 with Asm = 100 3 . 
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